Written by Aunoy Poddar July 9th, 2022
current_file <- rstudioapi::getActiveDocumentContext()$path
output_file <- stringr::str_replace(current_file, '.Rmd', '.R')
knitr::purl(current_file, output = output_file)
file.edit(output_file)
library(Seurat)
Warning message:
In fun(libname, pkgname) : couldn't connect to display ":0"
library(tictoc)
library(ggplot2)
library(patchwork)
library(pheatmap)
library(RColorBrewer)
library(tidyverse)
library(gridExtra)
library(png)
library(cowplot)
library(magick)
data_dir = '/home/aunoy/st/arc_profiling/st_analysis/hand_annotated_data/rethresholded'
meta_dir = '/home/aunoy/st/arc_profiling/st_analysis/hand_annotated_data/overlay'
output_dir_plot = '/home/aunoy/st/arc_profiling/st_analysis/results/plots'
output_dir_tbls = '/home/aunoy/st/arc_profiling/st_analysis/results/tables'
df_408 = data.frame()
for (file_name in list.files(data_dir)){
print(file_name)
if(grepl('164', file_name)){
next
}
#if(grepl('408_TC', file_name) | grepl('408_vMS', file_name)){
# next
#}
df_to_append <- read.table(file.path(data_dir, file_name), sep = ',', header = TRUE)
while(length(ind <- which(df_to_append$Image.Name == "")) > 0){
df_to_append$Image.Name[ind] <- df_to_append$Image.Name[ind -1]
}
colnames(df_to_append) <- toupper(colnames(df_to_append))
df_to_append <- df_to_append %>%
mutate(area = strsplit(file_name, '.csv')[[1]])
## Add relative_XY_position
if(!is_empty(df_408)){
df_to_append <- df_to_append %>%
dplyr::select(colnames(df_408))
}
df_408 <- rbind(df_408, df_to_append)
}
[1] "164_CC.csv"
[1] "164_MS_CC.csv"
[1] "164_MS_TC.csv"
[1] "164_TC.csv"
[1] "408_CC.csv"
[1] "408_dMS_TC.csv"
[1] "408_MS_CC.csv"
[1] "408_TC.csv"
[1] "408_vMS_TC.csv"
df_408$IMAGE.NAME = unlist(lapply(df_408$IMAGE.NAME, gsub, pattern='_Cluster', replacement=''))
df_408$IMAGE.NAME = unlist(lapply(df_408$IMAGE.NAME, gsub, pattern='[*]', replacement=''))
df_408$IMAGE.NAME = unlist(lapply(df_408$IMAGE.NAME, gsub, pattern='X', replacement=''))
df_408$IMAGE.NAME = unlist(lapply(df_408$IMAGE.NAME, gsub, pattern='L2_', replacement='L2-'))
df_408$IMAGE.NAME = unlist(lapply(df_408$IMAGE.NAME, gsub, pattern='-L2', replacement='_L2'))
df_408$IMAGE.NAME = unlist(lapply(df_408$IMAGE.NAME, gsub, pattern='Tc_12', replacement='TC_12'))
## Missing
df_408 = df_408[df_408$IMAGE.NAME != 'Layer1', ]
df_408 = df_408[df_408$IMAGE.NAME != 'TC_1', ]
df_408 = df_408[df_408$IMAGE.NAME != 'TC_18', ]
df_408 = df_408[df_408$IMAGE.NAME != 'TC_19', ]
#df_408$IMAGE.NAME = toupper(df_408$IMAGE.NAME)
unique(df_408$IMAGE.NAME)
[1] "CC_Cortical1" "CC_Cortical2" "CC_L2-1" "CC_L2-2" "CC_L2-3" "TC_2"
[7] "TC_3" "TC_4" "TC_5" "TC_6" "TC_7" "TC_8"
[13] "TC_9" "TC_10" "CC_4" "CC_5" "CC_6" "CC_7"
[19] "CC_8" "CC_9" "CC_10" "CC_11" "CC_12" "TC_16"
[25] "TC_17" "TC_20" "TC_11" "TC_12" "TC_13" "TC_14"
[31] "TC_15"
image_names = unique(df_408$IMAGE.NAME)
# Preset these variables to negative values so I can easily check if they were updated later
df_408$X = -1
df_408$Y = -1
# set some normalization variables
## This is the size of the image when the pixel values are taken from top left down
IMAGE_SIZE = 1024
## This is the size of an image in the global coordinate space
IMAGE_LEN = 20
# Load the dataframe with global and relative coordinates
img_cords = read.table(file.path(meta_dir, '408_pixel_coordinates.csv'), sep = ',', header = TRUE)
images = list.files(meta_dir)
for(image_name in image_names){
if(grepl('164', image_name)){
next
}
split_names = strsplit(image_name, '_')
cortex = toupper(split_names[[1]][1])
number = split_names[[1]][2]
number_csv = paste0('_', number, '.csv')
filename = images[grepl(cortex, images) & grepl(number_csv, images) & grepl('408', images)]
coordinates = read.table(file.path(meta_dir, filename), sep = ',', header = TRUE)
## checked already that lists are equal, missing 1, 18, 19 for now, layer 1 and others
if(cortex == 'CC'){
#print(paste('cc', filename, image_name))
x_adj = 0
y_adj = 0
} else if(as.numeric(number) <= 10){
#print(paste('tc<11', filename, image_name))
x_adj = 180#img_cords[img_cords$Name == 'G_TC_1', 'x']
y_adj = 210 #img_cords[img_cords$Name == 'G_TC_1', 'y']
}else{
#print(paste('tc>=11', filename, image_name))
x_adj = 380#img_cords[img_cords$Name == 'G_TC_11', 'x']
y_adj = 410 #img_cords[img_cords$Name == 'G_TC_11', 'y']
}
## so this is a little tricky, so need to get it right
## Remember, it is the top right that the coordinate is coming from, but
## the bottom right is the new coordinate space.
## so first when we get the original coordinate space, to set to relative
## of bottom would be the same X, but 1024 - Y
## push out the coordinates for better visualization
x_repelled <- (512 - coordinates$X_Coordinate_In_pixels)
df_408[df_408$IMAGE.NAME == image_name, 'X'] = (x_repelled /
IMAGE_SIZE * IMAGE_LEN) +
img_cords[img_cords$Name == image_name, 'x'] + x_adj
df_408[df_408$IMAGE.NAME == image_name, 'Y'] = ((1024-coordinates$Y_Coordinate_In_pixels) /
IMAGE_SIZE * IMAGE_LEN) +
img_cords[img_cords$Name == image_name, 'y'] + y_adj
}
jy_408 = df_408 %>%
dplyr::select(-c(area, IMAGE.NAME, X, Y)) %>%
t() %>%
CreateSeuratObject()
jy_408 <- NormalizeData(jy_408, scale.factor = 1e5) ###
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
normed = GetAssayData(jy_408, slot = 'data')
normed[normed < 3] = 0
jy_408 <- SetAssayData(jy_408, slot = 'data', normed)
xycords = df_408 %>% select(c('X', 'Y')) %>% as.matrix()
colnames(xycords) <- c('pixel_1', 'pixel_2')
jy_408[["XY"]] <- CreateDimReducObject(embeddings = xycords, key = "pixel_", assay = DefaultAssay(jy_408))
jy_408$gad1_true = normed['GAD1',] != 0 & normed['SATB1',] == 0
Error in intI(i, n = x@Dim[1], dn[[1]], give.dn = FALSE) :
invalid character indexing
nkx21 <- normed['NKX2.1',] != 0 & normed['LHX6',] == 0 & jy_408$gad1_true
lhx6 <- normed['NKX2.1',] == 0 & normed['LHX6',] != 0 & jy_408$gad1_true
mge_both <- normed['NKX2.1',] != 0 & normed['LHX6',] != 0 & jy_408$gad1_true
mge_neither <- normed['NKX2.1',] == 0 & normed['LHX6',] == 0 & jy_408$gad1_true
non_interneuron <- !jy_408$gad1_true
jy_408$mge_lineage = 'error'
jy_408$mge_lineage[nkx21] = 'NKX2.1'
jy_408$mge_lineage[lhx6] = 'LHX6'
jy_408$mge_lineage[mge_both] = 'NKX2.1 & LHX6'
jy_408$mge_lineage[mge_neither] = 'Non-MGE'
jy_408$mge_lineage[non_interneuron] = 'Non-IN'
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
cols = c('grey', 'purple'), reduction = "XY", pt.size = 0.2, group.by = 'gad1_true', order = which(jy_408$gad1_true))
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 0.1, split.by = 'mge_lineage',
cells.highlight = list(nkx2.1 = which(nkx21), lhx6 = which(lhx6), both = which(mge_both)),
cols.highlight = c('orange1','hotpink1', 'black'),
order = which(jy_408$gad1_true)) + scale_y_reverse() + scale_x_reverse() + NoAxes()
sp8 <- normed['SP8',] != 0 & normed['COUPTF2',] == 0 & jy_408$gad1_true
couptf2 <- normed['SP8',] == 0 & normed['COUPTF2',] != 0 & jy_408$gad1_true
cge_lge <- normed['SP8',] != 0 & normed['COUPTF2',] != 0 & jy_408$gad1_true
neither_cge_lge <- normed['SP8',] == 0 & normed['COUPTF2',] == 0 & jy_408$gad1_true
non_interneuron <- !jy_408$gad1_true
jy_408$cge_lineage = 'error'
jy_408$cge_lineage[sp8] = 'SP8'
jy_408$cge_lineage[couptf2] = 'COUPTF2'
jy_408$cge_lineage[cge_lge] = 'SP8 & COUPTF2'
jy_408$cge_lineage[neither_cge_lge] = 'Non-CGE/LGE'
jy_408$cge_lineage[non_interneuron] = 'Non-IN'
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 1, split.by = 'cge_lineage',
cells.highlight = list(sp8 = which(sp8), couptf2 = which(couptf2), both = which(cge_lge)),
cols.highlight = c('red','blue', 'purple'),
order = which(jy_408$gad1_true)) + scale_y_reverse() + scale_x_reverse() + NoAxes()
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('greenyellow', 'lightgoldenrodyellow', 'orange1', 'aquamarine4', 'tomato', 'dodgerblue3', 'violetred3'),
cols = c('greenyellow', 'lightgoldenrodyellow', 'aquamarine4', 'orange1', 'dodgerblue3', 'tomato', 'violetred3'),
reduction = "XY", order = c(6, 4, 3, 2, 0, 1, 5), pt.size = 1) + scale_y_reverse() + scale_x_reverse() #group.by = 'mge_lineage',
#cells.highlight = list(sp8 = which(sp8), couptf2 = which(couptf2), both = which(cge_lge)),
#cols.highlight = c('red','blue', 'purple'),
#order = which(jy_408$gad1_true)) + scale_y_reverse() + scale_x_reverse()
cluster_by_loc = as.data.frame(as.matrix(cbind(rownames(df_408), jy_408$seurat_clusters, df_408$X, df_408$Y)))
colnames(cluster_by_loc) = c('cellnum', 'cluster', 'X', 'Y')
cluster_by_loc %>%
#filter(cluster == 3) %>%
ggplot(aes(x = X, color = cluster)) + geom_histogram(position = 'identity', bins = 50, binwidth = 1000, stat = 'count')
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 1, #group.by = 'mge_lineage',
cells.highlight = list(cluster0 = which(jy_408$seurat_clusters == 0))) + scale_y_reverse() + scale_x_reverse()#,
#cols.highlight = c('red','blue', 'purple'),
#order = which(jy_408$gad1_true)) + scale_y_reverse() + scale_x_reverse()
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 1, #group.by = 'mge_lineage',
cells.highlight = list(cluster1 = which(jy_408$seurat_clusters == 1))) + scale_y_reverse() + scale_x_reverse()#,
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 1, #group.by = 'mge_lineage',
cells.highlight = list(cluster2 = which(jy_408$seurat_clusters == 2))) + scale_y_reverse() + scale_x_reverse()#,
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 1, #group.by = 'mge_lineage',
cells.highlight = list(cluster3 = which(jy_408$seurat_clusters == 3))) + scale_y_reverse() + scale_x_reverse()#,
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 1, #group.by = 'mge_lineage',
cells.highlight = list(cluster4 = which(jy_408$seurat_clusters == 4))) + scale_y_reverse() + scale_x_reverse()#,
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 1, #group.by = 'mge_lineage',
cells.highlight = list(cluster5 = which(jy_408$seurat_clusters == 5))) + scale_y_reverse() + scale_x_reverse()#,
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 1, #group.by = 'mge_lineage',
cells.highlight = list(cluster6 = which(jy_408$seurat_clusters == 6))) + scale_y_reverse() + scale_x_reverse()#,
FeaturePlot(jy_408, features = c('NKX2.1', 'LHX6'),
reduction = "XY", pt.size = 1, order = TRUE, split.by = 'area', by.col = TRUE) + scale_y_reverse() + scale_x_reverse() + NoAxes() + NoLegend()
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
cols = c('greenyellow', 'lightgoldenrodyellow', 'orange1', 'aquamarine4', 'tomato', 'dodgerblue3', 'violetred3'),
reduction = "XY", pt.size = 1, split.by = 'seurat_clusters') + scale_y_reverse() + scale_x_reverse() + NoAxes() + NoLegend()
jy_408.markers <- FindAllMarkers(jy_408, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
Calculating cluster 0
| | 0 % ~calculating
|+++++++++ | 17% ~00s
|+++++++++++++++++ | 33% ~00s
|+++++++++++++++++++++++++ | 50% ~00s
|++++++++++++++++++++++++++++++++++ | 67% ~00s
|++++++++++++++++++++++++++++++++++++++++++ | 83% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
Calculating cluster 1
| | 0 % ~calculating
|+++++ | 9 % ~00s
|++++++++++ | 18% ~00s
|++++++++++++++ | 27% ~00s
|+++++++++++++++++++ | 36% ~00s
|+++++++++++++++++++++++ | 45% ~00s
|++++++++++++++++++++++++++++ | 55% ~00s
|++++++++++++++++++++++++++++++++ | 64% ~00s
|+++++++++++++++++++++++++++++++++++++ | 73% ~00s
|+++++++++++++++++++++++++++++++++++++++++ | 82% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++ | 91% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
Calculating cluster 2
| | 0 % ~calculating
|+++++ | 10% ~00s
|++++++++++ | 20% ~00s
|+++++++++++++++ | 30% ~00s
|++++++++++++++++++++ | 40% ~00s
|+++++++++++++++++++++++++ | 50% ~00s
|++++++++++++++++++++++++++++++ | 60% ~00s
|+++++++++++++++++++++++++++++++++++ | 70% ~00s
|++++++++++++++++++++++++++++++++++++++++ | 80% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 90% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
Calculating cluster 3
| | 0 % ~calculating
|++++++ | 11% ~00s
|++++++++++++ | 22% ~00s
|+++++++++++++++++ | 33% ~00s
|+++++++++++++++++++++++ | 44% ~00s
|++++++++++++++++++++++++++++ | 56% ~00s
|++++++++++++++++++++++++++++++++++ | 67% ~00s
|+++++++++++++++++++++++++++++++++++++++ | 78% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 89% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
Calculating cluster 4
| | 0 % ~calculating
|++++ | 8 % ~00s
|++++++++ | 15% ~00s
|++++++++++++ | 23% ~00s
|++++++++++++++++ | 31% ~00s
|++++++++++++++++++++ | 38% ~00s
|++++++++++++++++++++++++ | 46% ~00s
|+++++++++++++++++++++++++++ | 54% ~00s
|+++++++++++++++++++++++++++++++ | 62% ~00s
|+++++++++++++++++++++++++++++++++++ | 69% ~00s
|+++++++++++++++++++++++++++++++++++++++ | 77% ~00s
|+++++++++++++++++++++++++++++++++++++++++++ | 85% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 92% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
Calculating cluster 5
| | 0 % ~calculating
|+++++ | 9 % ~00s
|++++++++++ | 18% ~00s
|++++++++++++++ | 27% ~00s
|+++++++++++++++++++ | 36% ~00s
|+++++++++++++++++++++++ | 45% ~00s
|++++++++++++++++++++++++++++ | 55% ~00s
|++++++++++++++++++++++++++++++++ | 64% ~00s
|+++++++++++++++++++++++++++++++++++++ | 73% ~00s
|+++++++++++++++++++++++++++++++++++++++++ | 82% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++ | 91% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
Calculating cluster 6
| | 0 % ~calculating
|++++++ | 11% ~00s
|++++++++++++ | 22% ~00s
|+++++++++++++++++ | 33% ~00s
|+++++++++++++++++++++++ | 44% ~00s
|++++++++++++++++++++++++++++ | 56% ~00s
|++++++++++++++++++++++++++++++++++ | 67% ~00s
|+++++++++++++++++++++++++++++++++++++++ | 78% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++ | 89% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
jy_408.markers %>%
group_by(cluster) %>%
slice_max(n = 32, order_by = avg_log2FC)
jy_408_IN <- jy_408[, jy_408$gad1_true]
jy_408_IN <- FindVariableFeatures(jy_408_IN, selection.method = "vst")
all.genes <- rownames(jy_408_IN)
jy_408_IN <- ScaleData(jy_408_IN, features = all.genes)
jy_408_IN <- RunPCA(jy_408_IN, approx = FALSE)
jy_408_IN <- FindNeighbors(jy_408_IN, dims = 1:30)
jy_408_IN <- FindClusters(jy_408_IN, resolution = 0.8)
jy_408_IN <- RunUMAP(jy_408_IN, dims = 1:30)
DimPlot(jy_408_IN, reduction = "umap", group.by = 'seurat_clusters')
jy_408_IN[["XY"]] <- CreateDimReducObject(embeddings = xycords[jy_408$gad1_true, ], key = "pixel_", assay = DefaultAssay(jy_408_IN))
DimPlot(jy_408_IN, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 1) + scale_y_reverse() + scale_x_reverse()
jy_408_IN.markers <- FindAllMarkers(jy_408_IN, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
jy_408_IN.markers %>%
group_by(cluster) %>%
slice_max(n = 32, order_by = avg_log2FC)
jy_408_IN <- jy_408[, jy_408$gad1_true]
jy_408_IN <- FindVariableFeatures(jy_408_IN, selection.method = "vst")
all.genes <- rownames(jy_408_IN)
jy_408_IN <- ScaleData(jy_408_IN, features = all.genes)
jy_408_IN <- RunPCA(jy_408_IN, approx = FALSE)
jy_408_IN <- FindNeighbors(jy_408_IN, dims = 1:30)
jy_408_IN <- FindClusters(jy_408_IN, resolution = 0.8)
jy_408_IN <- RunUMAP(jy_408_IN, dims = 1:30)
DimPlot(jy_408_IN, reduction = "umap", group.by = 'seurat_clusters')
jy_408_IN[["XY"]] <- CreateDimReducObject(embeddings = xycords[jy_408$gad1_true, ], key = "pixel_", assay = DefaultAssay(jy_408_IN))
DimPlot(jy_408_IN, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 1) + scale_y_reverse() + scale_x_reverse()
jy_408_IN.markers <- FindAllMarkers(jy_408_IN, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
jy_408_IN.markers %>%
group_by(cluster) %>%
slice_max(n = 32, order_by = avg_log2FC)
image_names = unique(df_408$IMAGE.NAME)
# Preset these variables to negative values so I can easily check if they were updated later
df_408$X = -1
df_408$Y = -1
# set some normalization variables
## This is the size of the image when the pixel values are taken from top left down
IMAGE_SIZE = 1024
## This is the size of an image in the global coordinate space
IMAGE_LEN = 20
# Load the dataframe with global and relative coordinates
img_cords = read.table(file.path(meta_dir, '408_pixel_coordinates.csv'), sep = ',', header = TRUE)
images = list.files(meta_dir)
for(image_name in image_names){
split_names = strsplit(image_name, '_')
cortex = toupper(split_names[[1]][1])
number = split_names[[1]][2]
number_csv = paste0('_', number, '.csv')
filename = images[grepl(cortex, images) & grepl(number_csv, images)]
coordinates = read.table(file.path(meta_dir, filename), sep = ',', header = TRUE)
## checked already that lists are equal, missing 1, 18, 19 for now, layer 1 and others
if(cortex == 'CC'){
#print(paste('cc', filename, image_name))
x_adj = 20
y_adj = 188
} else if(as.numeric(number) <= 10){
#print(paste('tc<11', filename, image_name))
x_adj = 410#img_cords[img_cords$Name == 'G_TC_1', 'x']
y_adj = 470 #img_cords[img_cords$Name == 'G_TC_1', 'y']
}else{
#print(paste('tc>=11', filename, image_name))
x_adj = 590#img_cords[img_cords$Name == 'G_TC_11', 'x']
y_adj = 790 #img_cords[img_cords$Name == 'G_TC_11', 'y']
}
## so this is a little tricky, so need to get it right
## Remember, it is the top right that the coordinate is coming from, but
## the bottom right is the new coordinate space.
## so first when we get the original coordinate space, to set to relative
## of bottom would be the same X, but 1024 - Y
## push out the coordinates for better visualization
#x_repelled <- (512 - coordinates$X_Coordinate_In_pixels)
df_408[df_408$IMAGE.NAME == image_name, 'X'] = (coordinates$X_Coordinate_In_pixels /
IMAGE_SIZE * IMAGE_LEN) +
img_cords[img_cords$Name == image_name, 'x']/1.5 + x_adj
df_408[df_408$IMAGE.NAME == image_name, 'Y'] = ((1024-coordinates$Y_Coordinate_In_pixels) /
IMAGE_SIZE * IMAGE_LEN) +
img_cords[img_cords$Name == image_name, 'y'] + y_adj
}
xycords = df_408 %>% select(c('X', 'Y')) %>% as.matrix()
colnames(xycords) <- c('pixel_1', 'pixel_2')
jy_408[["XY"]] <- CreateDimReducObject(embeddings = xycords, key = "pixel_", assay = DefaultAssay(jy_408))
#https://stackoverflow.com/questions/9917049/inserting-an-image-to-ggplot2
theme_set(theme_cowplot())
bad_colors <- DimPlot(jy_408, cols = c('greenyellow', 'lightgoldenrodyellow', 'aquamarine4', 'orange1', 'dodgerblue3', 'tomato', 'violetred3'), reduction = "XY", pt.size = 0.01, order = c(6, 4, 3, 2, 0, 1, 5)) + scale_y_reverse() + scale_x_reverse() + xlim(852, 0) + ylim(1242, 0) + NoLegend() + coord_fixed() + NoAxes()
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
#xorig = -852
#yorig = -1242
ggdraw() +
draw_image('~/st/arc_profiling/st_analysis/hand_annotated_data/images/408_slice_noboxes_nocolor.png',
x = 0, y = 0) +
draw_plot(bad_colors)
DimPlot(jy_408, reduction = "XY", pt.size = 0.1)+ scale_y_reverse() + scale_x_reverse() + xlim(852, 0) + ylim(1242, 0) # + NoLegend() + NoAxes()
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
jy_408_sp[["XY"]] <- CreateDimReducObject(embeddings = xycords, key = "pixel_", assay = DefaultAssay(jy_408_sp))
reln <- normed['VLDLR',] != 0 & normed['LRP8',] == 0 & jy_408$gad1_true
lrp8 <- normed['VLDLR',] == 0 & normed['LRP8',] != 0 & jy_408$gad1_true
relnboth <- normed['VLDLR',] != 0 & normed['LRP8',] != 0 & jy_408$gad1_true
relnneither <- normed['VLDLR',] == 0 & normed['LRP8',] == 0 & jy_408$gad1_true
non_interneuron <- !jy_408$gad1_true
jy_408$reln_signaling = 'error'
jy_408$reln_signaling[reln] = 'VLDLR'
jy_408$reln_signaling[lrp8] = 'LRP8'
jy_408$reln_signaling[relnboth] = 'VLDLR & LRP8'
jy_408$reln_signaling[relnneither] = 'Neither'
jy_408$reln_signaling[non_interneuron] = 'Non-IN'
## just define sets of cells that I want to plot
cingulate_MS = grepl('CC', df_408$IMAGE.NAME)
dorsal_TC_MS = df_408$IMAGE.NAME %in% outer('TC_', 2:10, FUN=paste0)
ventral_TC_MS = !cingulate_MS & !dorsal_TC_MS
FeaturePlot(jy_408, cells = which(ventral_TC_MS), reduction = "XY", features = 'TBR1') + scale_y_reverse() + scale_x_reverse()
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
FeaturePlot(jy_408, cells = which(ventral_TC_MS), reduction = "XY", features = 'LRP8') + scale_y_reverse() + scale_x_reverse()
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
FeaturePlot(jy_408, cells = which(dorsal_TC_MS), reduction = "XY", features = c('CXCR4', 'CALB2'), order = TRUE) + scale_y_reverse() + scale_x_reverse()
FeaturePlot(jy_408, cells = which(ventral_TC_MS), reduction = "XY", features = c('CXCR4', 'CALB2'), order = TRUE) + scale_y_reverse() + scale_x_reverse()
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
FeaturePlot(jy_408, cells = which(cingulate_MS), reduction = "XY", features = c('NKX2.1'), order = TRUE) + scale_y_reverse() + scale_x_reverse()
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.
FeaturePlot(jy_408, cells = which(ventral_TC_MS), reduction = "XY", features = c('CXCR4', 'CALB2'), order = TRUE) +
stat_density_2d(aes(fill = ..density..), geom = "raster", contour = FALSE) +
scale_fill_distiller(palette=4, direction=-1)
normed_data = t(as.matrix(GetAssayData(jy_408, slot = 'data')))
norm_bools = normed_data > 0
density_df = as.data.frame(cbind(xycords, norm_bools, jy_408$seurat_clusters))
colnames(density_df) = c('X', 'Y', colnames(norm_bools), 'cluster')
density_df$fov = 'error'
density_df$fov[dorsal_TC_MS] = 'dMS_TC'
density_df$fov[ventral_TC_MS] = 'vMS_TC'
density_df$fov[cingulate_MS] = 'MS_CC'
density_df %>%
ggplot(aes(x = X, y = Y)) +
geom_bin2d(bins = 20) +
scale_fill_continuous(type = "viridis") +
theme_bw()
DimPlot(jy_408, #cells = grepl('CC', df_408$area),
#cols = c('grey', 'purple', 'blue', 'red', 'hotpink1'),
reduction = "XY", pt.size = 0.1, split.by = 'reln_signaling',
cells.highlight = list(vldlr = which(reln), lrp8 = which(lrp8), both = which(relnboth)),
cols.highlight = c('orange1','hotpink1', 'black'),
order = which(jy_408$gad1_true)) + scale_y_reverse() + scale_x_reverse()
goi_primary = c('VIP', 'SST', 'CXCR4')
goi = c('VIP', 'SST', 'CXCR4', 'CXCR7', 'CXCL14', 'CXCL12', 'LRP8', 'VLDLR', 'RELN')
fxn_genes = as.data.frame(normed_data[, goi])
fxn_genes$cge = jy_408$cge_lineage
fxn_genes$fov = density_df$fov
fxn_genes %>%
filter(fov == 'dMS_TC') %>%
pivot_longer(!c(cge, fov), names_to = "gene", values_to = "expr") %>%
ggplot(aes(x = gene, y = expr, fill = cge)) + geom_boxplot()# geom_bar(position="dodge", stat="identity")
fxn_genes %>%
filter(fov == 'vMS_TC') %>%
pivot_longer(!c(cge, fov), names_to = "gene", values_to = "expr") %>%
ggplot(aes(x = gene, y = expr, fill = cge)) + geom_boxplot()# geom_bar(position="dodge", stat="identity")
cfov = 'MS_CC'
fxn_genes %>%
filter(fov == cfov) %>%
pivot_longer(!c(cge, fov), names_to = "gene", values_to = "expr") %>%
ggplot(aes(x = gene, y = expr, fill = cge)) + geom_bar(position="dodge", stat="identity") +
scale_fill_manual(values=c("green","gray","gray","blue", "red")) + ggtitle(cfov)
cfov = 'vMS_TC'
fxn_genes %>%
filter(fov == cfov) %>%
pivot_longer(!c(cge, fov), names_to = "gene", values_to = "expr") %>%
ggplot(aes(x = gene, y = expr, fill = cge)) + geom_bar(position="dodge", stat="identity") +
scale_fill_manual(values=c("green","gray","gray","blue", "red")) + ggtitle(cfov)
cfov = 'dMS_TC'
fxn_genes %>%
filter(fov == cfov) %>%
pivot_longer(!c(cge, fov), names_to = "gene", values_to = "expr") %>%
ggplot(aes(x = gene, y = expr, fill = cge)) + geom_bar(position="dodge", stat="identity") +
scale_fill_manual(values=c("green","gray","gray","blue", "red")) + ggtitle(cfov)
fxn_genes = as.data.frame(norm_bools[, goi])
fxn_genes$cge = jy_408$cge_lineage
fxn_genes$fov = density_df$fov
fxn_genes %>%
filter(fov == 'dMS_TC') %>%
pivot_longer(!c(cge, fov), names_to = "gene", values_to = "expr") %>%
ggplot(aes(x = gene, y = as.numeric(expr), fill = cge)) + geom_bar(position="dodge", stat="summary", fun = "mean") +
xlab("pct") + ggtitle('dMS_TC')
fxn_genes = as.data.frame(norm_bools[, goi])
fxn_genes$mge = jy_408$mge_lineage
fxn_genes$fov = density_df$fov
cfov = 'vMS_TC'
fxn_genes %>%
filter(fov == cfov) %>%
select(-fov) %>%
pivot_longer(!mge, names_to = "gene", values_to = "expr") %>%
data_summary(varname="expr", groupnames=c("gene", "mge")) %>%
ggplot(aes(x=gene, y=expr, fill=mge)) + geom_bar(stat="identity", color="black", position=position_dodge()) + geom_errorbar(aes(ymin=expr-sem, ymax=expr+sem), width=.2, position=position_dodge(.9)) +
ylab("pct") + ggtitle(cfov) + scale_fill_manual(values=c("cadetblue1","brown2","darkslateblue","grey", "khaki1")) +
theme_classic()
#c("cadetblue1","brown2","darkslateblue","grey", "khaki1") for mge
#http://www.sthda.com/english/wiki/ggplot2-error-bars-quick-start-guide-r-software-and-data-visualization
library(plotrix)
Attaching package: ‘plotrix’
The following object is masked from ‘package:fields’:
color.scale
data_summary <- function(data, varname, groupnames){
require(plyr)
summary_func <- function(x, col){
c(mean = mean(x[[col]], na.rm=TRUE),
sem = std.error(x[[col]], na.rm=TRUE))
}
data_sum<-ddply(data, groupnames, .fun=summary_func,
varname)
data_sum <- rename(data_sum, c("mean" = varname))
return(data_sum)
}
df2 <- fxn_genes %>%
select(-fov) %>%
pivot_longer(!cge, names_to = "gene", values_to = "expr") %>%
data_summary(varname="expr",
groupnames=c("gene", "cge"))
ggplot(df2, aes(x=gene, y=expr, fill=cge)) +
geom_bar(stat="identity", color="black", position=position_dodge()) +
geom_errorbar(aes(ymin=expr-sem, ymax=expr+sem), width=.2,
position=position_dodge(.9))
ggplot(df2, aes(x=gene, y=expr, fill=cge)) +
geom_bar(stat="identity", color="black", position=position_dodge()) +
geom_errorbar(aes(ymin=expr-sem, ymax=expr+sem), width=.2,
position=position_dodge(.9))
fxn_genes %>%
select(-c(fov, mge)) %>%
pivot_longer(!cge, names_to = "gene", values_to = "expr") %>%
data_summary(varname="expr", groupnames=c("gene", "cge")) %>%
ggplot(aes(x=gene, y=expr, fill=cge)) + geom_bar(stat="identity", color="black", position=position_dodge()) + geom_errorbar(aes(ymin=expr-sem, ymax=expr+sem), width=.2, position=position_dodge(.9)) +
ylab("pct") + scale_fill_manual(values=c("darkolivegreen1","orange2","rosybrown4","gray96", "grey")) +
theme_classic()
print(206/388)
[1] 0.5309278
unique(df_408$IMAGE.NAME)
[1] "CC_Cortical1" "CC_Cortical2" "CC_L2-1" "CC_L2-2" "CC_L2-3" "TC_2"
[7] "TC_3" "TC_4" "TC_5" "TC_6" "TC_7" "TC_8"
[13] "TC_9" "TC_10" "CC_4" "CC_5" "CC_6" "CC_7"
[19] "CC_8" "CC_9" "CC_10" "CC_11" "CC_12" "TC_16"
[25] "TC_17" "TC_20" "TC_11" "TC_12" "TC_13" "TC_14"
[31] "TC_15"
unique(df_408$IMAGE.NAME)
[1] "CC_Cortical1" "CC_Cortical2" "CC_L2-1" "CC_L2-2" "CC_L2-3" "TC_2"
[7] "TC_3" "TC_4" "TC_5" "TC_6" "TC_7" "TC_8"
[13] "TC_9" "TC_10" "CC_4" "CC_5" "CC_6" "CC_7"
[19] "CC_8" "CC_9" "CC_10" "CC_11" "CC_12" "TC_16"
[25] "TC_17" "TC_20" "TC_11" "TC_12" "TC_13" "TC_14"
[31] "TC_15"
orderd = c('TC_20', 'TC_17', 'TC_16', 'TC_15', 'TC_14', 'TC_13', 'TC_12', 'TC_11',
'TC_10', 'TC_9', 'TC_8', 'TC_7', 'TC_6', 'TC_4', 'TC_4', 'TC_3', 'TC_2',
'CC_4', 'CC_5', 'CC_6', 'CC_6', 'CC_8', 'CC_9', 'CC_10', 'CC_11', 'CC_12',
'CC_L2-1', 'CC_L2-2', 'CC_L2-3', 'CC_Cortical1', 'CC_Cortical2')
x_horz = 1:length(images_ordered) * 25
y_horz = rep(0, length(images_ordered))
horz_embedding = data.frame()
df_408$X_horz = -1
df_408$Y_horz = -1
images = list.files(meta_dir)
for(i in 1:length(images_ordered)){
image_name = images_ordered[i]
split_names = strsplit(image_name, '_')
cortex = toupper(split_names[[1]][1])
number = split_names[[1]][2]
number_csv = paste0('_', number, '.csv')
filename = images[grepl(cortex, images) & grepl(number_csv, images)]
coordinates = read.table(file.path(meta_dir, filename), sep = ',', header = TRUE)
## checked already that lists are equal, missing 1, 18, 19 for now, layer 1 and others
## so this is a little tricky, so need to get it right
## Remember, it is the top right that the coordinate is coming from, but
## the bottom right is the new coordinate space.
## so first when we get the original coordinate space, to set to relative
## of bottom would be the same X, but 1024 - Y
## push out the coordinates for better visualization
#x_repelled <- (512 - coordinates$X_Coordinate_In_pixels)
df_408[df_408$IMAGE.NAME == image_name, 'X_horz'] = (coordinates$X_Coordinate_In_pixels /
IMAGE_SIZE * IMAGE_LEN) + y_horz[i]
df_408[df_408$IMAGE.NAME == image_name, 'Y_horz'] = ((1024-coordinates$Y_Coordinate_In_pixels) /
IMAGE_SIZE * IMAGE_LEN) + x_horz[i]
}
hcoords = df_408 %>% dplyr::select(c('X_horz', 'Y_horz')) %>% as.matrix()
colnames(hcoords) <- c('pixel_1', 'pixel_2')
jy_408[["H"]] <- CreateDimReducObject(embeddings = hcoords, key = "pixel_", assay = DefaultAssay(jy_408))
Warning: Cannot add objects with duplicate keys (offending key: pixel_) setting key to original value 'h_'
jy_408$satb1_true = !jy_408$gad1_true
DimPlot(jy_408, cols = c('grey', 'purple'), reduction = "H", pt.size = 0.25, group.by = 'gad1_true', order = which(jy_408$gad1_true)) + coord_fixed(ratio = 0.5) + NoAxes() + NoLegend()
DimPlot(jy_408, reduction = "H", pt.size = 0.25, group.by = 'mge_lineage', order = "NKX2.1 & LHX6", cells.highlight = which(jy_408$mge_lineage == 'NKX2.1 & LHX6'), cols.highlight = 'red') + coord_fixed(ratio = 1) + NoLegend() + NoAxes()
theme_Publication <- function(base_size=14, base_family="helvetica") {
library(grid)
library(ggthemes)
(theme_foundation(base_size=base_size, base_family=base_family)
+ theme(plot.title = element_text(face = "bold",
size = rel(1.2), hjust = 0.5),
text = element_text(),
panel.background = element_rect(colour = NA),
plot.background = element_rect(colour = NA),
panel.border = element_rect(colour = NA),
axis.title = element_text(face = "bold",size = rel(1)),
axis.title.y = element_text(angle=90,vjust =2),
axis.title.x = element_text(vjust = -0.2),
axis.text = element_text(),
axis.line = element_line(colour="black"),
axis.ticks = element_line(),
panel.grid.major = element_line(colour="#f0f0f0"),
panel.grid.minor = element_blank(),
legend.key = element_rect(colour = NA),
legend.position = "bottom",
legend.direction = "horizontal",
legend.key.size= unit(0.2, "cm"),
legend.margin = unit(0, "cm"),
legend.title = element_text(face="italic"),
plot.margin=unit(c(10,5,5,5),"mm"),
strip.background=element_rect(colour="#f0f0f0",fill="#f0f0f0"),
strip.text = element_text(face="bold")
))
}
scale_fill_Publication <- function(...){
library(scales)
discrete_scale("fill","Publication",manual_pal(values = c("#386cb0","#fdb462","#7fc97f","#ef3b2c","#662506","#a6cee3","#fb9a99","#984ea3","#ffff33")), ...)
}
scale_colour_Publication <- function(...){
library(scales)
discrete_scale("colour","Publication",manual_pal(values = c("#386cb0","#fdb462","#7fc97f","#ef3b2c","#662506","#a6cee3","#fb9a99","#984ea3","#ffff33")), ...)
}
FeaturePlot(jy_408, pt.size=0.1, reduction = "H", features = 'LRP8') + coord_fixed(ratio = 0.5) + theme_Publication() + NoLegend() + NoAxes()
Warning: `legend.margin` must be specified using `margin()`. For the old behavior use legend.spacing